snATACseq data processing

This dataset contains nuclei of CD45+CD11b+F4/80+CD115+ FACS sorted monocyte enriched cells from mouse E14 fetal liver and adult bone marrow. The input files can be downloaded from GSE292830.

The mouse brain Signac vignette was followed for the data processing: https://stuartlab.org/signac/articles/mouse_brain_vignette

suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(Signac))
suppressPackageStartupMessages(library(Seurat, lib.loc="/home/daliya/Apps/Seurat.v4"))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(cowplot))
suppressPackageStartupMessages(library(clustree))
suppressPackageStartupMessages(library(plotly))
suppressPackageStartupMessages(library(presto))
suppressPackageStartupMessages(library(scater))
suppressPackageStartupMessages(library(EnsDb.Mmusculus.v79))

Set the working directory

knitr::opts_knit$set(root.dir = "~/Documents/Kia/Jonathan/JBA1-2_ATAC/results/")

path<-"/media/daliya/SSD2/0-0-Exp4885-JBA/"
path.output<-"~/Documents/Kia/Jonathan/JBA1-2_ATAC/results/"
sample.names<-c("JBA1","JBA2")
sample.description=c("BM","Fetal liver")
mm10.blacklist = rtracklayer:: import(gzfile(  "/media/daliya/SSD2/ATAC_blacklist/mm10.blacklist.bed.gz")) ### Blacklist downloaded from https://github.com/Boyle-Lab/Blacklist?tab=readme-ov-file

Pre-processing

When pre-processing chromatin data, Signac uses information from two related input files, both of which can be created using CellRanger:

• Peak/Cell matrix. This is analogous to the gene expression count matrix used to analyze single-cell RNA-seq. However, instead of genes, each row of the matrix represents a region of the genome (a peak), that is predicted to represent a region of open chromatin. Each value in the matrix represents the number of Tn5 integration sites for each single barcode (i.e. a cell) that map within each peak. You can find more detail on the 10X Website.
• Fragment file. This represents a full list of all unique fragments across all single cells. It is a substantially larger file, is slower to work with, and is stored on-disk (instead of in memory). However, the advantage of retaining this file is that it contains all fragments associated with each single cell, as opposed to only fragments that map to peaks. More information about the fragment file can be found on the 10x Genomics website or on the sinto website.

We start by creating a Seurat object using the peak/cell matrix and cell metadata generated by cellranger-atac, and store the path to the fragment file on disk in the Seurat object:

seur=list()
for (i in 1:length(sample.names)){
chrassay<- CreateChromatinAssay(
  counts =  Read10X_h5(paste0(path,sample.names[i],"/filtered_peak_bc_matrix.h5")),
  sep = c(":", "-"),
  genome = "mm10",
  fragments = paste0(path,sample.names[i],'/fragments.tsv.gz'),
  min.cells = 1,
  verbose =F
)
metadata <- read.csv(paste0(path,sample.names[i],"/singlecell.csv"),  header = TRUE,row.names = 1)
seur[[i]] <- CreateSeuratObject(
  counts = chrassay,
  assay = 'peaks',
  project = 'ATAC',
  meta.data = metadata
)
seur[[i]]$sample=sample.names[i]
seur[[i]]$origin=sample.description[i]
}
rm(chrassay)
rm(metadata)

We can also add gene annotations to the brain object for the mouse genome. This will allow downstream functions to pull the gene annotation information directly from the object.

annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Mmusculus.v79, verbose=F)

# change to UCSC style since the data was mapped to hg19
seqlevels(annotations) <- paste0('chr', seqlevels(annotations))
genome(annotations) <- "mm10"

# add the gene information to the object

for (i in 1:length(sample.names)){
Annotation(seur[[i]]) <- annotations
}

Computing QC Metrics

for (i in 1:length(sample.names)){
seur[[i]] <- NucleosomeSignal(object = seur[[i]], verbose=F)
}

• Nucleosome banding pattern We can look at the fragment length periodicity for all the cells, and group by cells with high or low nucleosomal signal strength. Cells which are outliers for the mononucleosomal/ nucleosome-free ratio have different banding patterns. The remaining cells exhibit a pattern that is typical for a successful ATAC-seq experiment.

for (i in 1:length(sample.names)){
seur[[i]]$nucleosome_group <- ifelse(seur[[i]]$nucleosome_signal > 4, 'NS > 4', 'NS < 4')
print(FragmentHistogram(object = seur[[i]], group.by = 'nucleosome_group', region = 'chr1-1-10000000'))
}

• TSS enrichment score The enrichment of Tn5 integration events at transcriptional start sites (TSSs) can also be an important quality control metric to assess the targeting of Tn5 in ATAC-seq experiments. The ENCODE consortium defined a TSS enrichment score as the number of Tn5 integration site around the TSS normalized to the number of Tn5 integration sites in flanking regions. See the ENCODE documentation for more information about the TSS enrichment score (https://www.encodeproject.org/data-standards/terms/).

for (i in 1:length(sample.names)){
seur[[i]] <- TSSEnrichment(seur[[i]], verbose=F)
}

• Fraction of fragments in peaks Represents the fraction of all fragments that fall within ATAC-seq peaks. Cells with low values (i.e. <15-20%) often represent low-quality cells or technical artifacts that should be removed. Note that this value can be sensitive to the set of peaks used.

for (i in 1:length(sample.names)){
seur[[i]] $pct_reads_in_peaks <- seur[[i]] $peak_region_fragments / seur[[i]] $passed_filters * 100
}

• Total number of fragments in peaks: “peak_region_fragments” metric, calculated by CellRanger A measure of cellular sequencing depth / complexity. Cells with very few reads may need to be excluded due to low sequencing depth. Cells with extremely high levels may represent doublets, nuclei clumps, or other artefacts.

• Ratio reads in genomic blacklist regions The ENCODE project has provided a list of blacklist regions, representing reads which are often associated with artefactual signal. Cells with a high proportion of reads mapping to these areas (compared to reads mapping to peaks) often represent technical artifacts and should be removed.

seur[[i]]$blacklist_ratio <- FractionCountsInRegion(
  object = seur[[i]], 
  assay = 'peaks',
  regions = mm10.blacklist
)

Fidning outlier nuclei

The exact QC thresholds for removing outliers need to be adjusted according to the dataset. Here, we use the isOutlier() function from the scuttle package. It defines an observation as an outlier if it is more than a specified number of median absolute deviations (MADs, default 3) from the median in the specified direction.

  j= 'nucleosome_signal'
for (i in 1:length(sample.names)){
  outliers= scuttle::isOutlier(seur[[i]]@meta.data[,j], nmads=6, type="higher", log=TRUE)
  attr(outliers, "thresholds") <- NULL
 seur[[i]]@meta.data[paste0(j,".outlier.higher")] <- outliers
  cat(sample.names[i],j,"outliers higher :",sum(seur[[i]]@meta.data[paste0(j,".outlier.higher")]),"\n")
}

JBA1 nucleosome_signal outliers higher : 218 
JBA2 nucleosome_signal outliers higher : 223 

for ( j in c('pct_reads_in_peaks','TSS.enrichment')){
 for (i in 1:length(sample.names)){
    outliers= scater::isOutlier(seur[[i]]@meta.data[,j], nmads=6, type="lower", log=TRUE)
  attr(outliers, "thresholds") <- NULL
 seur[[i]]@meta.data[paste0(j,".outlier.lower")] <- outliers
  cat(sample.names[i],j,"outliers lower :",sum(seur[[i]]@meta.data[paste0(j,".outlier.lower")]),"\n")
  }
}

JBA1 pct_reads_in_peaks outliers lower : 326 
JBA2 pct_reads_in_peaks outliers lower : 214 
JBA1 TSS.enrichment outliers lower : 19 
JBA2 TSS.enrichment outliers lower : 55 

j='peak_region_fragments'
for (i in 1:length(sample.names)){
    outliers= scater::isOutlier(seur[[i]]@meta.data[,j], nmads=2, type="lower", log=TRUE)
  attr(outliers, "thresholds") <- NULL
 seur[[i]]@meta.data[paste0(j,".outlier.lower")] <- outliers
  cat(sample.names[i],j,"outliers lower :",sum(seur[[i]]@meta.data[paste0(j,".outlier.lower")]),"\n")
}

JBA1 peak_region_fragments outliers lower : 859 
JBA2 peak_region_fragments outliers lower : 1412 

j="blacklist_ratio"
for (i in 1:length(sample.names)){
  outliers= scater::isOutlier(seur[[i]]@meta.data[,j], nmads=3, type="higher", log=TRUE)
  attr(outliers, "thresholds") <- NULL
 seur[[i]]@meta.data[paste0(j,".outlier.higher")] <- outliers
  cat(sample.names[i],j,"outliers higher :",sum(seur[[i]]@meta.data[paste0(j,".outlier.higher")]),"\n")
}

JBA1 blacklist_ratio outliers higher : 35 
JBA2 blacklist_ratio outliers higher : 315 

for ( j in c('peak_region_fragments','pct_reads_in_peaks', 'nucleosome_signal',
               'TSS.enrichment', "blacklist_ratio")){
  for (i in 1:length(sample.names)){
     hist(seur[[i]]@meta.data[,j],
        breaks = 100,xlab=j,
        main=paste0(j,": ",sample.names[i]))
    if (paste0(j,".outlier.higher") %in% colnames(seur[[i]]@meta.data) ){
      if(sum(seur[[i]]@meta.data[paste0(j,".outlier.higher")]) !=0 )
       abline(v = min(seur[[i]]@meta.data[,j][unlist(seur[[i]]@meta.data[paste0(j,".outlier.higher")])]), col = "red")
    }
   if (paste0(j,".outlier.lower") %in% colnames(seur[[i]]@meta.data) ){
    if(sum(seur[[i]]@meta.data[paste0(j,".outlier.lower")]) !=0 )
   abline(v = max(seur[[i]]@meta.data[,j][unlist(seur[[i]]@meta.data[paste0(j,".outlier.lower")])]), col = "red")
   }
  }
}

We remove cells that are outliers for these QC metrics.

for (i in 1:length(sample.names)){
seur[[i]]  <- seur[[i]] [,!(seur[[i]]$peak_region_fragments.outlier.lower | 
    seur[[i]]$pct_reads_in_peaks.outlier.lower | 
    seur[[i]]$nucleosome_signal.outlier.higher  | 
    seur[[i]]$TSS.enrichment.outlier.lower |
   seur[[i]]$blacklist_ratio.outlier.higher   )]

print(sample.names[i])
print(seur[[i]] )
}

[1] "JBA1"
An object of class Seurat 
178086 features across 18322 samples within 1 assay 
Active assay: peaks (178086 features, 0 variable features)
 2 layers present: counts, data
[1] "JBA2"
An object of class Seurat 
180414 features across 17338 samples within 1 assay 
Active assay: peaks (180414 features, 0 variable features)
 2 layers present: counts, data

seur<-merge(seur[[1]],seur[2:length(seur)])

Normalization and linear dimensional reduction

Normalization: Signac performs term frequency-inverse document frequency (TF-IDF) normalization. This is a two-step normalization procedure, that both normalizes across cells to correct for differences in cellular sequencing depth, and across peaks to give higher values to more rare peaks.

Feature selection: The low dynamic range of scATAC-seq data makes it challenging to perform variable feature selection, as we do for scRNA-seq. Instead, we can choose to use only the top n% of features (peaks) for dimensional reduction, or remove features present in less than n cells with the FindTopFeatures() function. Here we will use all features, though we have seen very similar results when using only a subset of features (try setting min.cutoff to ‘q75’ to use the top 25% all peaks), with faster runtimes. Features used for dimensional reduction are automatically set as VariableFeatures() for the Seurat object by this function.

Dimension reduction: We next run singular value decomposition (SVD) on the TD-IDF matrix, using the features (peaks) selected above. This returns a reduced dimension representation of the object (for users who are more familiar with scRNA-seq, you can think of this as analogous to the output of PCA).

The combined steps of TF-IDF followed by SVD are known as latent semantic indexing (LSI), and were first introduced for the analysis of scATAC-seq data by Cusanovich et al. 2015.

seur <- RunTFIDF(seur,verbose =F)
seur <- FindTopFeatures(seur, min.cutoff = 'q0',verbose =F)
seur <- RunSVD(object = seur,verbose =F)

The first LSI component often captures sequencing depth (technical variation) rather than biological variation. If this is the case, the component should be removed from downstream analysis. We can assess the correlation between each LSI component and sequencing depth using the DepthCor() function:

DepthCor(seur)

ElbowPlot(object = seur,ndims =50, reduction = "lsi")

Non-linear dimension reduction and clustering

Now that the cells are embedded in a low-dimensional space we can use methods commonly applied for the analysis of scRNA-seq data to perform graph-based clustering and non-linear dimension reduction for visualization. The functions RunUMAP(), FindNeighbors(), and FindClusters() all come from the Seurat package.

### LSI components selection for downstream analysis
dims.use<-30

seur <- RunUMAP(object = seur, reduction = 'lsi', dims = 2:dims.use,verbose =F)
seur <- FindNeighbors(seur,  reduction = 'lsi',dims = 2:dims.use, verbose=F, graph.name=paste0("ATAC_snn_LSI",dims.use))
for ( i in seq(0,2, 0.25))
  seur <- FindClusters(seur, resolution = i, algorithm = 3, graph.name=paste0("ATAC_snn_LSI",dims.use), verbose=F) 

clustree::clustree(seur, prefix = paste0("ATAC_snn_LSI",dims.use,"_res."))+
  ggtitle(paste("LSI =",dims.use))

plot<-list()
for ( res in c(0.25,0.5,0.75,1))
  plot[[as.character(res)]]<-DimPlot(seur,label=T,repel=T, group.by = paste0("ATAC_snn_LSI",dims.use,"_res.",res))+
                                      ggtitle(paste("LSI =",dims.use,"res=",res))
plot_grid(plotlist=plot)

DimPlot(seur,repel =T,label=T, group.by = "origin") 

Create a gene activity matrix

We can try to quantify the activity of each gene in the genome by assessing the chromatin accessibility associated with the gene, and create a new gene activity assay derived from the scATAC-seq data. Here we will use a simple approach of summing the fragments intersecting the gene body and promoter region

To create a gene activity matrix, we extract gene coordinates and extend them to include the 2 kb upstream region (as promoter accessibility is often correlated with gene expression). We then count the number of fragments for each cell that map to each of these regions, using the using the FeatureMatrix() function. These steps are automatically performed by the GeneActivity() function:

gene.activities <- GeneActivity(seur,verbose =F)

seur[['RNA']] <- CreateAssayObject(counts = gene.activities)
 # add the gene activity matrix to the Seurat object as a new assay and normalize it                               
seur <- NormalizeData(
  object = seur,
  assay = 'RNA',
  normalization.method = 'LogNormalize',
  scale.factor = median(seur$nCount_RNA),
  verbose =F
)
seur <- ScaleData(seur, assay = 'RNA',verbose =F)

DefaultAssay(seur) <- 'RNA'
res=0.5
Idents(seur)=  paste0("ATAC_snn_LSI",dims.use,"_res.",res)
Idents(seur)=  factor(Idents(seur),levels = 0:(length(unique(Idents(seur)))-1))
DEgenes=list()
for (i in levels(Idents(seur))){
 DEgenes[[i]]<-FindMarkers(seur, ident.1 = i,min.cells.group=2,pseudocount.use = 0.01, max.cells.per.ident = 1000)
 DEgenes[[i]]$cluster=i
 DEgenes[[i]]$pseudocount=0.01
 DEgenes[[i]]$max.cells.per.ident=1000
 DEgenes[[i]]$gene=row.names(DEgenes[[i]])
}

features.use=unlist(lapply(DEgenes, function(x) { head(x[x$avg_log2FC>0,]$gene,7)}))
length(features.use)

[1] 77

names(features.use)=NULL
features.use=features.use[!duplicated(features.use)]

DimPlot(seur,repel =T,label=T) 

DotPlot(seur, features = features.use)+RotatedAxis()+NoLegend()

Idents(seur) <- plyr::mapvalues(x = Idents(seur), from = 0:10, 
    to =c("FL 2", "BM 3","BM 1","FL 1","BM 4","BM 2",
          "FL 3", "FL 5", "Erythroid","FL 4","Lymphoid"))

new.order=c (  "FL 1", "FL 2" , "FL 3" ,  "FL 4","FL 5" ,  "BM 1"  , "BM 2", "BM 3"   ,  "BM 4" ,  "Erythroid" , "Lymphoid"       )
new.order[!new.order %in% levels( Idents(seur))]

character(0)

levels( Idents(seur))[!levels( Idents(seur)) %in%new.order]

character(0)

Idents(seur)=factor( Idents(seur), levels=new.order)
seur$annot=Idents(seur)

DimPlot(seur,repel =T,label=T) 

DotPlot(seur, features=c("Ptprc","Flt3","Tshz2","Igf2bp3" , "Gli2","Cd34","Prtn3","Kit","Mpo","Lpo","Elane" ,"Sorcs2","Pax5","Ror1","Tmtc1","C1qa","C1qb","Ms4a7","Apoe","Fcrls","Cdk8","Lars2","Plec","S100a10","F13a1","Ly6c2","Crip1" ,"Ccr2" ,"Lyz2", "Tgfbi","Mrc1","Ust","H2-Aa","Cd74","Atf3" ,"Klf13"   , "Grk5"  ,"Dock5"  , "Cyp2ab1" ,  "Ace","Ear2","Ank1","Gypa","Aqp1","Car2","Hbb-bt" ,"Ccr7","Cd79a","Ms4a1","Ebf1","Dntt","Mki67","Stmn1","nucleosome_signal","nucleosome_percentile","TSS.enrichment","TSS.percentile","pct_reads_in_peaks","blacklist_ratio"))+RotatedAxis()+theme(axis.title = element_blank())

meta.data<-seur@meta.data %>%
    dplyr::group_by(annot,origin) %>%
    dplyr::summarise(count=dplyr::n())%>% suppressMessages()%>%
    dplyr::group_by(annot) %>%
    dplyr::mutate(perc.per.group = (count / sum(count))*100)  
ggplot(meta.data, aes(x=annot, y=perc.per.group, fill=origin))+
    geom_bar(stat = "identity")+
    theme_classic()+
    theme(axis.title.x=element_blank(),
          axis.text = element_text(size=12),
          axis.title.y = element_text(size=12),
          legend.text=element_text(size=12),
          legend.title=element_text(size=12))+
    ylab("% cells per cluster")+RotatedAxis()

saveRDS(seur,paste0("Integrated.BM.and.FL_monocytes_snATACseq_JBA1-2.seurat.rds"))

sessionInfo()

R version 4.4.3 (2025-02-28)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=de_BE.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=de_BE.UTF-8   
 [6] LC_MESSAGES=en_US.UTF-8    LC_PAPER=de_BE.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=de_BE.UTF-8 LC_IDENTIFICATION=C       

time zone: Europe/Brussels
tzcode source: system (glibc)

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] EnsDb.Mmusculus.v79_2.99.0  ensembldb_2.30.0            AnnotationFilter_1.30.0     GenomicFeatures_1.58.0     
 [5] AnnotationDbi_1.68.0        scater_1.34.0               scuttle_1.16.0              SingleCellExperiment_1.28.1
 [9] SummarizedExperiment_1.36.0 Biobase_2.66.0              GenomicRanges_1.58.0        GenomeInfoDb_1.42.1        
[13] IRanges_2.40.0              S4Vectors_0.44.0            BiocGenerics_0.52.0         MatrixGenerics_1.18.0      
[17] matrixStats_1.4.1           presto_1.0.0                data.table_1.16.2           Rcpp_1.0.13-1              
[21] plotly_4.10.4               clustree_0.5.1              ggraph_2.2.1                cowplot_1.1.3              
[25] ggplot2_3.5.1               SeuratObject_5.0.2          Seurat_4.3.0                Signac_1.14.0              
[29] dplyr_1.1.4                

loaded via a namespace (and not attached):
  [1] RcppAnnoy_0.0.22         splines_4.4.3            later_1.4.1              BiocIO_1.16.0            bitops_1.0-9            
  [6] tibble_3.2.1             polyclip_1.10-7          XML_3.99-0.17            lifecycle_1.0.4          globals_0.16.3          
 [11] lattice_0.22-5           MASS_7.3-65              backports_1.5.0          magrittr_2.0.3           sass_0.4.9              
 [16] rmarkdown_2.29           jquerylib_0.1.4          yaml_2.3.10              httpuv_1.6.15            sctransform_0.4.1       
 [21] spam_2.11-0              sp_2.1-4                 spatstat.sparse_3.1-0    reticulate_1.40.0        DBI_1.2.3               
 [26] pbapply_1.7-2            RColorBrewer_1.1-3       abind_1.4-8              zlibbioc_1.52.0          Rtsne_0.17              
 [31] purrr_1.0.2              RCurl_1.98-1.16          tweenr_2.0.3             GenomeInfoDbData_1.2.13  ggrepel_0.9.6           
 [36] irlba_2.3.5.1            listenv_0.9.1            spatstat.utils_3.1-3     goftest_1.2-3            spatstat.random_3.3-3   
 [41] fitdistrplus_1.2-1       parallelly_1.40.0        leiden_0.4.3.1           codetools_0.2-20         DelayedArray_0.32.0     
 [46] RcppRoll_0.3.1           ggforce_0.4.2            tidyselect_1.2.1         UCSC.utils_1.2.0         farver_2.1.2            
 [51] ScaledMatrix_1.14.0      viridis_0.6.5            spatstat.explore_3.4-2   GenomicAlignments_1.42.0 jsonlite_1.8.9          
 [56] BiocNeighbors_2.0.1      tidygraph_1.3.1          progressr_0.15.1         ggridges_0.5.6           survival_3.8-3          
 [61] tools_4.4.3              ica_1.0-3                glue_1.8.0               gridExtra_2.3            SparseArray_1.6.0       
 [66] xfun_0.49                withr_3.0.2              fastmap_1.2.0            fansi_1.0.6              digest_0.6.37           
 [71] rsvd_1.0.5               R6_2.5.1                 mime_0.12                colorspace_2.1-1         scattermore_1.2         
 [76] tensor_1.5               RSQLite_2.3.8            spatstat.data_3.1-6      utf8_1.2.4               tidyr_1.3.1             
 [81] generics_0.1.3           rtracklayer_1.66.0       graphlayouts_1.2.1       httr_1.4.7               htmlwidgets_1.6.4       
 [86] S4Arrays_1.6.0           uwot_0.2.2               pkgconfig_2.0.3          gtable_0.3.6             blob_1.2.4              
 [91] lmtest_0.9-40            XVector_0.46.0           htmltools_0.5.8.1        dotCall64_1.2            ProtGenerics_1.38.0     
 [96] scales_1.3.0             png_0.1-8                spatstat.univar_3.1-2    knitr_1.49               rstudioapi_0.17.1       
[101] rjson_0.2.23             reshape2_1.4.4           checkmate_2.3.2          curl_6.0.1               nlme_3.1-167            
[106] zoo_1.8-12               cachem_1.1.0             stringr_1.5.1            KernSmooth_2.23-26       vipor_0.4.7             
[111] parallel_4.4.3           miniUI_0.1.1.1           restfulr_0.0.15          pillar_1.9.0             grid_4.4.3              
[116] vctrs_0.6.5              RANN_2.6.2               promises_1.3.2           BiocSingular_1.22.0      beachmat_2.22.0         
[121] xtable_1.8-4             cluster_2.1.8.1          beeswarm_0.4.0           evaluate_1.0.1           cli_3.6.3               
[126] compiler_4.4.3           Rsamtools_2.22.0         rlang_1.1.4              crayon_1.5.3             future.apply_1.11.3     
[131] labeling_0.4.3           plyr_1.8.9               ggbeeswarm_0.7.2         stringi_1.8.4            viridisLite_0.4.2       
[136] deldir_2.0-4             BiocParallel_1.40.0      munsell_0.5.1            Biostrings_2.74.0        lazyeval_0.2.2          
[141] spatstat.geom_3.3-6      Matrix_1.7-1             patchwork_1.3.0          bit64_4.5.2              future_1.34.0           
[146] KEGGREST_1.46.0          shiny_1.9.1              ROCR_1.0-11              igraph_2.1.2             memoise_2.0.1           
[151] bslib_0.8.0              fastmatch_1.1-4          bit_4.5.0.1